sf::st_read function also work? Why ?CLUES
Use the readRDS function.
library(sf)
iris_31 <- readRDS("../data/iris_31.rds")
# iris_31 <- st_read("../data/iris_31.rds")readRDS function.
plot(iris_31). What do you notice ?plot(iris_31)
# We notice that R performs 3 graphs: one graph per variable in the sf file.sf::st_geometry function? What solution do you propose then?# ?sf::st_geometry
# This function makes it possible to isolate the geography information of the sf object so that we put aside other variables (here CODE_IRIS, P14_POP and INSEE_COM).
plot(st_geometry(iris_31))
CLUES
Use the sf::st_crs function to guess the projection and sf::st_transform to change the projection.
#?st_crs
st_crs(iris_31)
par(mar = c(0,0,0,0), mfrow = c(1,2))
plot(st_geometry(iris_31))
iris_31_aeqd <- st_transform(iris_31, crs="+proj=aeqd +lat_0=90 +lon_0=0")
plot(st_geometry(iris_31_aeqd))
mat <- st_distance(x = iris_31[1:5,], y = iris_31[1:5,])
matUnits: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00 51153.341 56509.22 53137.481 13756.98
[2,] 51153.34 0.000 20957.21 3086.011 49828.91
[3,] 56509.22 20957.212 0.00 11077.630 61863.18
[4,] 53137.48 3086.011 11077.63 0.000 54345.46
[5,] 13756.98 49828.910 61863.18 54345.458 0.00
mat_aeqd <- st_distance(x = iris_31_aeqd[1:5,], y = iris_31_aeqd[1:5,])
mat_aeqdUnits: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00 57206.922 62306.32 59385.205 13798.28
[2,] 57206.92 0.000 20957.48 3204.855 55368.23
[3,] 62306.32 20957.481 0.00 11076.220 66744.65
[4,] 59385.20 3204.855 11076.22 0.000 59843.30
[5,] 13798.28 55368.225 66744.65 59843.302 0.00
identical(mat,mat_aeqd)[1] FALSE
# The two matrices are differentCLUES
Use the classic functions of the dplyr package: select, group_by et summarize. These functions also work with sf objects.
library(dplyr)
com_31 <- iris_31 %>%
select(INSEE_COM,P14_POP) %>%
group_by(INSEE_COM) %>%
summarize(P14_POP= sum(P14_POP)) %>%
st_cast("MULTIPOLYGON")
plot(st_geometry(com_31))
sir_31 <- readRDS("../data/sir_31.rds")
com_31 <- left_join(com_31,
sir_31 %>%
mutate(INSEE_COM=paste0(DEPET,COMET)) %>%
group_by(INSEE_COM) %>%
summarize(nb_of_rest= n()) %>%
st_set_geometry(NULL),
by=c("INSEE_COM"="INSEE_COM")
) %>%
mutate(nb_of_rest=ifelse(is.na(nb_of_rest),0,nb_of_rest))table_MAUP <- readRDS("../data/table_MAUP.rds") %>%
select(CODGEO,EPCI)
epci_31 <- com_31 %>%
left_join(table_MAUP,by=c("INSEE_COM"="CODGEO")) %>%
group_by(EPCI) %>%
summarize(P14_POP=sum(P14_POP),nb_of_rest= sum(nb_of_rest)) %>%
st_cast("MULTIPOLYGON")
plot(st_geometry(epci_31))
CLUES
Start by creating the layers EPCI_less100 and COM_more100 corresponding respectively to the EPCIs which contain less than 100 restaurants and to the municipalities which belong to an EPCI containing more than 100 restaurants. Then, use the do.call(rbind, list(EPCI_less100,COM_more100)) statement to merge them.
EPCI_less100 <- epci_31 %>% filter(nb_of_rest < 100) %>%
setNames(c("territory","P14_POP","nb_of_rest","geometry"))
COM_more100 <- left_join(com_31,table_MAUP,by=c("INSEE_COM"="CODGEO")) %>%
filter(!EPCI%in%EPCI_less100$territory) %>%
select(-EPCI) %>%
setNames(c("territory","P14_POP","nb_of_rest","geometry"))
par(mar = c(0,0,0,0), mfrow = c(1,2))
plot(st_geometry(epci_31))
plot(st_geometry(EPCI_less100), col="lightblue",add=TRUE)
plot(st_geometry(epci_31))
plot(st_geometry(COM_more100), col="lightblue",add=TRUE)
mix_31 <- do.call(rbind, list(EPCI_less100,COM_more100))
plot(st_geometry(mix_31), col="lightblue")


CLUES
The propSymbolsLayer function allows you to draw proportional circles.
library(cartography)
plot(st_geometry(mix_31), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
propSymbolsLayer(mix_31, var = "nb_of_rest", inches = 0.2)
getBreaks et carto.pal functions of the {cartography} package. For the creation of the “typo” variable, you can use the cut function and apply the parameters digit.lab = 2 and include.lowest = TRUE.
library(sf)
library(cartography)
fra <- st_read("../data/fra.shp", quiet = TRUE)
mix_31 <- readRDS("../data/mix_31.rds")
mix_31 <- mix_31 %>% mutate(nb_rest_10000inhab = 10000 * nb_of_rest / P14_POP)
bks <- getBreaks(v = mix_31$nb_rest_10000inhab, method = "quantile", nclass = 6)[-c(1:2)]
cols <- carto.pal("orange.pal", 4)
mix_31 <- mix_31 %>% mutate(typo = cut(nb_rest_10000inhab,breaks = bks, dig.lab = 2, include.lowest = TRUE))

With cartography :
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(mix_31), border="azure")
plot(st_geometry(fra), col="ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)],add=TRUE)
choroLayer(mix_31, var = "nb_rest_10000inhab", breaks = bks, col = cols, border = "grey80", lwd = 0.4,
legend.pos = "topleft", legend.title.txt = "Number of restaurants\nfor 10000 inhabitants",
add = TRUE)
propSymbolsLayer(mix_31, var="nb_of_rest", col=viridis::viridis(1,alpha=0.8),border=NA, legend.pos="left", legend.title.txt = "Number of restaurants", add = TRUE)
layoutLayer(title = "Restaurants", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)With ggplot2 :
library(ggplot2)
map_ggplot <- ggplot() +
geom_sf(data = fra, colour = "ivory3",
fill = "ivory") +
geom_sf(data = mix_31, aes(fill = typo), colour = "grey80") +
scale_fill_manual(name = "Number of restaurants\nfor 10000 inhabitants",
values = cols, na.value = "#303030")+
geom_sf(data = mix_31 %>% st_centroid(),
aes(size= nb_of_rest), color = viridis::viridis(1,alpha=0.8), show.legend = 'point')+
scale_size(name = "Number of restaurants",
breaks = c(1, 500, 2000),
range = c(0,18))+
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(mix_31)[c(1,3)],
ylim = st_bbox(mix_31)[c(2,4)]
) +
theme_minimal() +
theme(panel.background = element_rect(fill = "azure",color=NA)) +
labs(
title = "Restaurants",
caption = "Insee, 2018\nKim & Tim, 2018"
)
plot(map_ggplot)Pour Timothée, élements des fonctions fufun et fufun2 à modifier :
mapview function.
library(mapview)
library(sf)
library(cartography)
sir_31 <- readRDS("../data/sir_31.rds")
mapview(sir_31, map.types = "OpenStreetMap", col.regions = "#940000", label = paste(sir_31$L2_NORMALISEE, sir_31$NOMEN_LONG, sep = " - "), color = "white", legend = TRUE, layer.name = "Restaurants in SIRENE", homebutton = FALSE, lwd = 0.5)